Metabasin dynamics and local structure in supercooled water 
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We employ the Distance Matrix method to investigate metabasin dynamics in supercooled water. 
We find that the motion of the system consists in the exploration of a finite region of configuration 
space (enclosing several distinct local minima), named metabasin, followed by a sharp crossing to 
a different metabasin. The characteristic time between metabasin transitions is comparable to the 
structural relaxation time, suggesting that these transitions are relevant for the long time dynamics. 
The crossing between metabasins is accompanied by very rapid diffusional jumps of several groups 
of dynamically correlated particles. These particles form relatively compact clusters and act as 
cooperative relaxing units responsible for the density relaxation. We find that these mobile particles 
are often characterized by an average coordination larger than four, i.e. are located in regions where 
the tetrahedral hydrogen bond network is distorted. 

PACS numbers: 61.20.Ja; 61.20.Lc; 64.70.Pf 



I. INTRODUCTION 



The elucidation of the physical underpinnings of 
glassy dynamics represents a central open problem 
in condensed matter [lll3,IH. Within this realm, the 
discovery that glassy systems are dynamically het- 
erogeneous represented a major breakthrough and 
reinforced the interest in the Adam & Gibbs hetero- 
geneous scenario for glassy relaxation: the existence 
of cooperatively relaxing regions in the sample whose 
timescales and sizes grow as the temperature is de- 
creased in the supercooled regime [J d, H, 0, B Hi • 
A complimentary and powerful framework for this 
field has been provided by the so-called landscape 
paradigm d, [TO, [ll|, [l^, 113]: the dynamics of a 
glassy system can be followed as the exploration it 
performs of its Potential Energy Landscape (PEL), 
the surface generated by the potential energy of the 
system as a function of particle coordinates. The 
PEL can be partitioned into disjoint basins, where a 
basin is unambiguously defined as the set of points 
in configuration space connected to the same lo- 
cal minimum (called inherent structure, IS) via a 
steepest-descendant trajectory. Closely related IS 
are arranged in metabasins (MB). At low enough 
temperature the dynamics of the system can be de- 
composed in fast vibrations within IS and structural 
rearrangement transitions between different IS. In 
this terms, the dynamics of the system can be stud- 
ied following its dynamics between IS. If the tem- 
perature is very low the system doesn't have enough 



energy to easily overcome barriers between PEL val- 
leys, this resulting in a slowing down of the dynam- 
ics. The link between these two descriptions has 
been already studied by means of molecular dynam- 
ics simulations @, [H, [H, [l^ [l^. An interesting 
finding in such direction for binary Lennard- Jones 
systems suggested a compelling picture for glassy dy- 
namics 8]: the a- relaxation corresponds to a small 
number of fast crossings from one metabasin to a 
neighboring one involving the collective motion of 
a significant number of particles that form a rela- 
tively compact cluster, called "democratic cluster" 
or d-cluster. Since the influence of the structure 
on dynamics is only local in time [13], the long 
time dynamics could be described as a random walk 
on metabasins [13] triggered by d-clusters [l7| . 
Thus, it would be of great interest to learn whether 
these results hold valid for other glassy systems be- 
yond binary Lennard- Jones systems so as to deter- 
mine the generality of such a scenario. In this sense, 
the aim of the present work is to study the possi- 
ble metabasin structure of supercooled water and 
its connections to local structure and dynamics. 



II. METABASIN STRUCTURE IN 
SUPERCOOLED WATER 



Water is not only important from its condition of be- 
ing our most familiar liquid. Supercooled water ex- 



hibits an unusual behavior that has fostered intense 
experimental and theoretical work both in thermo- 
dynamics and dynamics in the past decades [H, 
[iqI . [20I [2l| . From a dynamical viewpoint, it has 
been shown computationally that supercooled wa- 
ter confirms the general picture of dynamical het- 
erogeneities '9', mils [23| , with mobile molecules ar- 
ranged in clusters usually linked by hydrogen bonds. 
The properties of its PEL have also been studied, in- 
cluding the characterization of transitions between 
ISs [lal- Within this context, it has also been 
demonstrated that the presence of molecules of high 
coordination (more than the usual four first neigh- 
bors demanded by its local, long-range disordered, 
tetrahedral structure) and of bifurcated hydrogen 
bonds (hydrogen bonds that link one hydrogen to 
two oxygens, instead of only one) enhance the local 
mobility [H [H, |2|. That regions of greater lo- 
cal density and the presence of bifurcated hydrogen 
bonds aid molecular mobility is not surprising since 
the possibility for higher coordination and hydrogen 
bonding can help lowering the barrier to restructure 
the hydrogen bond network. Indeed, for such water 
molecules, a new hydrogen bond can begin to be 
formed at the same time as the old one is being bro- 
ken, thus lowering the energetic cost [l^ . 

In this work we have studied a system of 216 water 
molecules interacting via the simple point charge ex- 
tended (SPC/E) model H in the (NVE) ensemble. 
The integration step is 1 ft second, and long range 
interactions are taken in account using the reaction 
field method. At temperatures within the super- 
cooled regime we investigated both the real dynam- 
ics (the instantaneous MD configurations) and the 
inherent dynamics, the dynamics of the quenched 
structures, which have been calculated performing 
periodical quenches of configurations obtained in the 
real trajectories using a standard conjugate gradient 
minimization algorithm with 10^^^ tolerance [23|. In 
order to study the possible metabasin structure of 
this system we recorded 250 configurations at fixed 
time intervals of 10 ps — at density 1 g/cm^ and 
T=210 K — . The total run time is larger than the 
a-relaxation. For these conditions the a-relaxation 
time has been estimated (as the time when the in- 
termediate incoherent scattering function, measured 
at the wave vector of the first peak in the structure 
factor, has decayed to 1/e) to be Tq = 769ps [28j , 
both for the real and the inherent dynamics. We 
shall present results for a run at T = 210 K but the 
same behavior was found for other trajectories at 
such temperature and at other close temperatures. 
With the recorded configurations we built the fol- 
lowing distance matrix @, [2^, DM: A'^{t',t") = 

:^Z]ilikj(^') - rj(^")l^) where ri{t) is the posi- 



tion of the oxygen of molecule i at time t (in th? 
real or the inherent configuration). Thus A'^{t',t") 
gives the system averaged squared displacement of 
a molecule in the time interval that starts at t' and 
ends at t" . In other words, this distance matrix con- 
tains the averaged squared distances between (the 
oxygens of) each recorded configuration and all the 
other ones. For this kind of study, we must investi- 
gate small systems, since for large systems different 
subsystems would produce interfering results [1, [l3| . 
Thus, we used N = 216 molecules (however, the 
same qualitative results would be obtained for a 
small subsystem immersed in a big system, thus 
ruling out the possibility for finite size effects, as 
found for binary Lennard- Jones systems Q). For 
this small-sized system it has been shown that no 
significant finite-size effects are present in the stud- 
ied temperature range 28]. . Fig. [T] shows typical 
results for a run at T = 210 K. The gray level of the 
squares in the DM depicts the distance between the 
corresponding configurations, the darker the shad- 
ing indicating the lower the distance between them. 
Fig. [T^ shows the real dynamics case while Fig. [T]d 
shows the inherent dynamics one. Additionally, in 
Fig. [TJ: we plot the hydrogen bonds Hamming dis- 
tance matrix, that is, we calculated the elements of 
this matrix as the number of hydrogen bond changes 
between the corresponding structures at times t' and 
t" (in the inherent dynamics) . For the identification 
of hydrogen bonds, HB, we employed a geometri- 
cal criterion consistent with previous mixed geomet- 
ric, energetic and dynamic strategies [30l|: two water 
molecules are linked by a HB if the 0-0 distance is 
lower than 0.35 nm (a distance that falls beyond the 
first peak of the 0-0 radial distribution function) 
and the 0-H...0 angle is lower than 60°. Neverthe- 
less, the vast majority of linear HBs found for the 
inherent structures presented 0-0 distances lower 
than 0.3 nm and very low angles (angles of 40° - 
60° occurred only for one of the bonds of bifurcated 
HBs, as we shall see later on). Direct inspection 
of these graphs shows that the results of the real 
and the inherent dynamics are practically identical 
(thus, from now on, we shall refer further results 
to the inherent dynamics, unless indication on the 
contrary). From the island structure of these matri- 
ces a clear MB structure of the landscape is evident. 
That is, islands are made up of closely related or 
similar configurations (low A^) which are separated 
from the configurations of other islands by large dis- 
tances. From the size of the islands found for many 
different independent runs at T = 210/^ wc can esti- 
mate the typical residence time in the MBs for this 
temperature to be close to the a-relaxation time, 
Tq,. Given the small system size we expect this to 
be a good estimate (however, this timescale clearly 



depends on system size, since for a large system dif- 
ferent subsystems would be undergoing MB transi- 
tion events at different times). We also studied the 
case of T = 220K and T = 230K (always for sys- 
tems of 216 molecules) for which Tq is respectively 
around 8 and 60 times lower than that for T = 210K 
and found that for both cases the typical MB resi- 
dence time is also close to the corresponding value 
of Tq,. Thus, Ta and the MB residence time seem to 
follow the same temperature dependence. The tran- 
sitions between the MBs (which, as can be learnt 
from Fig. [H last typically a few tens of picoseconds) 
are fast events compared to the times for the explo- 
ration of the MBs. Thus, consistent with the situa- 
tion for Lennard- Jones systems Q, the MB transi- 
tions are the events that trigger the a-relaxation in 
supercooled water. The hydrogen bonds Hamming 
DM showed in Fig. [TJ:, which is very similar to the 
DMs of Fig. [1^ and b, makes evident the fact that 
MB transitions imply significant restructuring of the 
hydrogen bonds (HB) network of supercooled water. 



III. METABASIN DYNAMICS AND 
MOLECULAR MOBILITY 



A still open question is if MB transitions are due 
to the large displacement of a few water molecules 
or to the result of more extensive collective rear- 
rangements, as that found in binary Lennard-Jones 
systems Q. To answer this point we performed a 
detailed analysis of molecular mobility, as shown in 
Fig. [21 The thick solid line of such figure represents 
the average squared displacement plot A^(0, t), that 
is, the first row of the DM of Fig. W^a.) Q . The thin 
solid line depicts the average squared displacement 
of the molecules (by monitoring the displacements of 
the oxygen atoms) in time interval 9— 40 ps, S'^{t, 6), 
which corresponds to A'^{t',t") measured along the 
diagonaH" =t'+e j3]. The bars of Fig.[l]display the 
function m(i, 9) Q, which represents the fraction of 
molecules that have moved more than a threshold of 
0.06 nm in time intervals [t,t + 9] with 6—10 ps (this 
threshold value is arbitrary but the reason for its 
choice will be given later on, and similar results arise 
for other values). From the thick solid line of Fig. [5] 
we can learn that MB transitions represent abrupt 
jumps (of 5-10 times the value of the squared dis- 
tances between contiguous configurations separated 
by 10 ps). In turn, the thin solid line makes clear the 
fact that MB transitions imply a significant enhance- 
ment of molecular mobility, while the bars relate this 
enhancement of global mobility to the fraction of 
molecules moving. Thus, in agreement with previous 



findings for binary Lennard-Jones systems @, M^ 
transitions imply the rearrangement of a substan- 
tial amount of the molecules of the system. That is, 
our results clearly demonstrate that also for super- 
cooled water the a-relaxation time corresponds to a 
small number of crossings from one metabasin to a 
neighboring one, each crossing being very rapid and 
involving the collective motion of a great number 
of particles. From direct inspection of the function 
■m{t, 9) in Fig. [5] we can see that in the MB tran- 
sitions 35 — 45% of the particles move more than 
0.06nrn in lOps, a time span slightly larger than 1% 
Ta (for comparison, an integration of the self part 
of the van Hove function, ATTr'^Gs{r,d), from 0.06 to 
infinity gives less than 0.15, cf. Fig. d]). The events 
that constitute the local exploration of the MB rep- 
resent the plateau regions of the thick solid line of 
Fig. [5] and, while relevant to the /?- relaxation, do not 
contribute to the structural relaxation of the system. 
Thus, these results demonstrate the prevailing role 
of metabasin transitions for the long time dynamics. 
Fig. [3]sohd line depicts n//B(0, t),the (number) dif- 
ference in HBs between a configuration at time t and 
the first configuration, while the bars indicate the 
function nHsit, t + 9), the number of HB that have 
changed between successive recorded configurations 
(separated by a time interval of length 9). The simi- 
larity between this figure and Fig.[2]imphes that the 
HB can also be used to signal MB transitions since 
such events entail a major rearrangement of the hy- 
drogen bond network of the system. In Fig. 3] we 
study the range of mobility of the molecules associ- 
ated with MB transitions. The plot of Fig.[3^ corre- 
sponds to the Van Hove function 47rr^Gs(r, t,t + 6) 
Q evaluated at chosen MB transition events (within 
[t,t + 9] time intervals, 6=10 ps) of the run under 
study. In other words. Fig. displays the distri- 
bution of mobility for the water molecules (by look- 
ing at the oxygens) at such events. The solid line 
depicts the distribution function of oxygen displace- 
ments averaged over the whole run, that is, the van 
Hove function ATrr'^Gs{r,6) (which gives the prob- 
ability for an oxygen to be located at a distance r 
from the origin after time 6) . Fig. is the analog of 
Fig. Hi but for the function Airr'^Gsir, t,t + 6) eval- 
uated at selected time intervals within MBs. From 
these plots we can learn that molecular motion is 
greatly enhanced at the MB transitions since the 
displacement distribution is not shifted compared to 
the van Hove function in Fig. HJd and it is signifi- 
cantly shifted to the right in Fig. [4^. Additionally, 
we can see that there is no clear preference for dis- 
placements of any given length. On average, the 
distribution functions for the MB transitions exceed 
the van Hove function at lengths r around 0.025 nm 
- 0.03 nm. Thus, we could select mobile molecules 



as those whose displacement is larger than such cut 
off. If we had used such value for the bars of Fig. [2] 
we should have obtained the same qualitative shape. 
However, we used a larger value (0.06 nm) since this 
will be more practical for the following section where 
we shall display and study the clusters of mobile 
particles at MB transitions. We also calculated the 
same quantities shown in Fig. 2^ for the real dy- 
namics — as opposed to the IS dynamics — where a 
similar behavior is found but with both kind of func- 
tions displaced to the right due to the thermal mo- 
tion. These data are not reported in this manuscript. 
While both approaches produce clear results, these 
are neater for the inherent dynamics. 



IV. DEMOCRATIC CLUSTERS AND 
LOCAL WATER STRUCTURE 



Having established the main role of large scale 
molecular rearrangements to the MB dynamics and 
thus, on the long time dynamics of supercooled wa- 
ter, the aim of this section is to study the nature 
of such molecular motions. Thus, the spatial distri- 
bution of mobile molecules and its relation to the 
local structure of the system will be explored. We 
shall see that, consistent with the situation for bi- 
nary Lennard- Jones systems @, relatively compact 
clusters of molecules conforming cooperative relax- 
ing units emerge as responsible for the structural 
relaxation of the system. Moreover, we shall make 
evident the connection of these clusters with local 
structural "defects" such as highly coordinated wa- 
ter molecules and bifurcated hydrogen bonds, effects 
that ha ve p reviously been regarded as mobility pro- 
moters O [13, nil- Fig. [S] shows a three-dimensional 
plot of the mobile particles (mobility greater than 
0.06 nm) for one of the time intervals of Fig. [H 
namely that which goes from t=1320 ps to ^=1330 
ps, but similar results were found for other MB tran- 
sitions. We show each mobile particle (oxygens and 
hydrogens) indicating the position at t=1320 ps and 
attaching a vector that shows its displacement to the 
position occupied at time i=1330 ps. A clear cluster 
arrangement of the mobile water molecules is evident 
from such picture. In this graph we can see that the 
mobile molecules are located close to the borders of 
the simulation box and an empty space occurs at 
center (the molecules in such region are not mobile) . 
These kind of compact clusters resemble the situa- 
tion in Lennard- Jones systems where were termed 
as "democratic clusters" or d-clusters @. These d- 
clusters represent potential candidates for the co- 
operatively relaxing regions proposed long ago by 



Adam & Gibbs whose verification has remained elu- 
sive till now (a d-cluster would represents a event 
when a subsystem, that reached a state that permit- 
ted a rearrangement, has performed such transition, 
that is, it would constitute an active subsystem at 
such time). However, we note that no definite con- 
nection between the d-clusters and the CRRs can be 
done without an exhaustive study of the dependence 
of the size of the d-clusters with temperature and its 
relation to the structural relaxation time. The mo- 
tions of many of the molecules of the d-clusters seem 
to be highly coordinated, with clear coherent flows of 
groups of molecules (they move in an organized fash- 
ion, that is, carrying a common course, as can be ver- 
ified by looking at the vectors attached to each atom 
in the 3D plot of the d-cluster of Fig. [5]) . To establish 
possible links between local structure and dynam- 
ics, we studied the local environment of the mobile 
molecules that comprise the d-clusters. An inter- 
esting connection between structure and dynamics 
has been found previously for supercooled water [isj 
and has been shown to be related to transitions be- 
tween contiguous ISs of the PEL which entail the 
movement of a group of a few water molecules: mo- 
bility is facilitated by the presence of defects like 
highly coordinated molecules (molecules with more 
than four first neighbors) and molecules involved in 
bifurcated hydrogen bonds (molecules which are hy- 
drogen bonded to two other water molecules via the 
same hydrogen). These regions of higher local den- 
sity aid mobility by lowering the cost of hydrogen 
bond reformulation since the breakage of the old 
hydrogen bond implied in the motion of the mo- 
bile molecule is somehow counterbalanced by the si- 
multaneous formation of the new one (a fact which 
would be otherwise anti- intuitive) [isj . Regarding 
bifurcated hydrogen bonds, we found that almost 
in all cases one of the bonds is almost linear (very 
low 0-H...0 angle) and short, while the other one 
is a bit longer and with a less favorable 0-H...0 an- 
gle between 40° and 60° (thus, this bond is clearly 
very weak). We note that minimized configurations 
present bifurcated HBs, a fact that speaks of cer- 
tain stability of such "species" which can thus be re- 
garded as intermediates rather than transition states 
for the interchange of hydrogen bonds. We shall 
show here results for the inherent dynamics, but the 
situation is similar for the real dynamics (albeit the 
identification of the HBs is blurred by thermal mo- 
tion resulting in a larger number of bifurcated HBs) . 
Fig. 6a shows the water molecules that posses more 
than four neighbors within a distance of 0.30 nm to- 
gether with such neighbors, for the structure at time 
i=1320 ps (the first structure of the MB transition 
interval shown in Fig. [5]). In turn. Fig. shows 
the bifurcated HBs present at the same configura- 



tion. We can easily learn that both highly coordi- 
nated molecules and bifurcated HBs are located at 
the same regions of the mobile molecules of Fig. [5] 
while the central region of the box, where no mo- 
bile particles occur, is also empty. That is, a clear 
link exists between the regions of high coordination 
and bifurcated bonds in a structure previous to a 
MB transition event, and the region where the im- 
minent d-cluster occurs, thus relating local structure 
to dynamics. For convenience, we shall term the wa- 
ter molecules that display coordination higher than 
four, the corresponding neighbors of such highly co- 
ordinated molecules and the molecules involved in 
bifurcated HBs as "defects" (for example, all the 
molecules displayed in Figs [5^ and b). We shall also 
define the function n(t, 9) as the number of coinci- 
dences or matches between the defects of the con- 
figuration at time t and the mobile molecules (with 
mobility greater than 0.06 nm) for the time inter- 
val [t, t + 9]. The outcome of this quantity is shown 
in Fig. [71 We can easily verify that this function is 
related to the MB transitions by noting the coinci- 
dences between the high bars of Fig. [7] and that of 
Fig. [2] (that indicates the previously defined func- 
tion 'm{t, 9) which monitors the number of mobile 
molecules in the corresponding time interval). Thus, 
the function n(t, 9) can also be used to signal MB 
transitions. 



ration of a finite region of configuration space iden- 
tified as metabasins followed by a sharp crossing to 
a different metabasin. The characteristic time be- 
tween metabasin transitions is comparable to the 
structural relaxation time. These transitions are 
therefore relevant for the long time dynamics. The 
crossing between metabasins is accompanied by very 
rapid diffusional jumps of several groups of dynam- 
ically correlated particles. These particles form rel- 
atively compact clusters and act as cooperative re- 
laxing units responsible for the density relaxation. 
These d-clusters represent potential candidates for 
the cooperatively relaxing regions proposed long ago 
by Adam & Gibbs. Finally, we related the mobility 
of the particles in the d-clusters to a studied quan- 
tity, namely the average local coordination number: 
mobile particles are located in regions where the 
tetrahedral hydrogen bond network is distorted. 
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V. CONCLUSIONS 



In this paper we have shown that the dynamics of 
water in the supercooled region consist in the explo- 
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FIG. 1: (a) and (b): Distance matrix A'^{t',t") re- 
spectively for real and inherent dynamics, (c): Hy- 
drogen bond Hamming distance matrix of the system 
for T = 210 K. The gray level correspond to values of 
A^(t',t") that arc given to the right of the figures (units 
are for (a) and (b) and hydrogen bond changes for 
(c)) 
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FIG. 2: Left scale: Average squared displacement 
A^(0,t) (thick solid line) and 5^{t,9) (thin solid line) 
for the trajectory given in Fig. The value of 9 is 
40 ps for this last curve. For both curves the units are 
nm^. Vertical bars (right scale): The function m{t,6) 
which gives the fraction of oxygen atoms that moved 
more than the threshold value rth ~ 0.06 nm in the time 
interval [t, t + 9], using 9 = 10 ps. 




FIG. 3: Solid line (left scale): Hydrogen bonds (num- 
ber) changes between a configuration at time t and the 
first configuration. Vertical bars (right scale): Hydrogen 
bonds (number) that have changed between successive 
recorded configurations (separated by a time interval of 
length 10 ps). 
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FIG. 4: The distribution function Gsir,t,t + 6) for dif- 
ferent values of t (curves with symbols). The value of 
9 is 10 ps. The bold curve is ATir'^Gsij-,9), the self part 
of the van Hove function, a) Values of t for which the 
system is about to leave a MB. b) Values of t in which 
the system is inside a MB. 



FIG. 5: Configuration snapsiiot of mobile molecules oc- 
curring in the MB-MB transition t = 1320 ps to f = 1330 
ps. The spheres (light and daxk for the hydrogen and 
oxygen atoms, respectively ) give the location of the 
atoms of the water molecule before the rearrangement 
and the arrows point to the position of the oxygen atom 
after the transition (to help visualize the particle rnov- 
ments, these arrows are twice the length of the corre- 
sponding displacement vectors. 



FIG. 6: Configuration snapshot of (a): molecules hav- 
ing more than four first neighbors within a distance of 
0.30 nm together with such neighbors (the central water 
molecule is fully depicted with its corresponding hydro- 
gen atoms, while only the oxygens of the first neighbors 
are included); (b): molecules with bifurcated hydrogen 
bonds, for the structure at time t = 1320 ps. Only 
the oxygens of the corresponding molecules arc indicated 
together with the hydrogens involved in the bifurcated 
bonds. Light and dark spheres are for the hydrogen and 
oxygen atoms, respectively. 
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FIG. 7: The function n(t, 6) shows the number of coiuci- 
dences or matches between the defects at time t and the 
mobile molecules (with mobility greater than 0.06 nm) 
in the time interval [t, t + 9] , with ^ = 10 ps. 



